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Numerical  calculation  of  gravity-capillary  interfacial  waves 
of  finite  amplitude 

Jean-Marc  Vanden-Broeck 

Department  of  Mathematics.  Stanford  University,  Stanford,  California  94305 
(Received  26  February  1980;  accepted  20  June  1980) 

Gravity -capillary  progressive  interfacial  waves  on  the  interface  between  two  semi-infinite  fluids  of  different 
densities  are  considered.  An  integro-differential  equation  for  the  unknown  shape  of  the  interface  is  derived. 

By  introducing  a  mesh  and  finite  difference,  this  equation  is  converted  into  a  finite  set  of  nonlinear  algebraic 
equations.  These  equations  are  solved  by  Newton's  method.  In  the  limiting  case  of  pure  gravity  waves,  the 
results  obtained  are  in  good  agreement  with  previous  calculations.  Two  continuous  families  of  capillary- 
gravity  waves  are  studied.  A  generalization  of  Wilton's  ripples  for  interfacial  waves  is  presented. 


I.  INTRODUCTION 

Over  the  past  decade  important  progress  had  been 
achived  in  the  calculation  of  free  surface  waves.  Sch¬ 
wartz1  extended  Stokes’  series  for  pure  gravity  waves 
to  high  order  by  computer  and  then  recast  these  series 
as  Pad£  approximants.  High  accuracy  solutions  were 
obtained  in  that  way.  Since  then,  this  technique  has  been 
applied  successfully  to  different  kinds  of  surface  waves 
such  as  solitary  waves  and  gravity -capillary  waves.2"5 

On  the  other  hand,  a  number  of  investigators9*9  have 
obtained  high  accuracy  solutions  by  using  direct  num¬ 
erical  approaches  based  on  an  integro-differential  equa¬ 
tion  formulation.  These  calculations  confirm  the  valid¬ 
ity  of  the  use  of  the  Padd  approximants  as  applied  to 
surface  waves.  However,  for  very  steep  waves  the 
numerical  approach  turns  out  to  be  more  efficient  than 
the  Padd  approximant  technique. 

Very  little  work  has  been  done  on  the  computation  of 
nonlinear  interfacial  waves.  The  main  results  are  con¬ 
tained  in  a  paper  by  Holyer10  who  used  Padd  approxi¬ 
mants  to  compute  pure  gravity  interfacial  waves. 

Inthe  present  paper,  we  compute  gravity- capillary  in¬ 
terfacial  waves  by  a  direct  numerical  scheme.  In  Sec. 
n,  the  problem  is  formulated  as  a  nonlinear  integro- 
differential  equation  for  the  unknown  shape  of  the  in¬ 
terface.  In  Sec.  Ill,  this  integral  equation  is  solved  by 
Newton’s  iterations.  The  mathematical  formulation 
and  the  numerical  method  are  similar  in  philosophy,  if 
not  in  detail,  to  the  procedure  used  by  Schwartz  and 
Vanden-Broeck7  and  Vanden-Broeck  and  Schwartz8  to 
compute  surface  waves.  The  results  are  presented 
in  Sec.  IV.  In  the  particular  case  of  pure  gravity  waves 
our  results  agree  with  those  obtained  by  Holyer.10 
Thus ,  the  validity  of  the  use  of  the  Padd  table  as  applied 
to  interfacial  waves  is  confirmed.  Many  different  fam¬ 
ilies  of  gravity-capillary  interfacial  waves  exist.  In 
Sec.  IVB  we  shall  study  two  of  them,  A  direct  gener¬ 
alization  of  Wilton’s  ripples  is  presented. 

II.  MATHEMATICAL  FORMULATION 

We  consider  two-dimensional  progressive  waves  of 
wavelength  x  and  phase  velocity  c,  propagating  under  the 
combined  effect  of  gravity  g  and  surface  tension  7 
along  the  interface  between  two  fluids.  We  measure 
lengths  in  units  of  X  and  velocities  in  units  of  c ,  so 


that  all  variables  become  dimensionless.  We  also 
choose  a  frame  of  reference  in  which  the  flows  are 
steady. 

The  fluids  are  assumed  to  be  incompressible  and 
irrotational.  Thus,  we  define  stream  functions  t%  and 
ip2  and  potential  functions  <px  and  <p2  for  the  lower  and 
upper  fluids,  respectively.  Without  loss  of  generality 
we  choose  =  i}2  =  0  on  the  interface  and  <t>x  =  <f>2  =  0 
at  one  crest.  Next,  we  introduce  rectangular  coordin¬ 
ates  (x,y)  with  the  x  axis  parallel  to  the  velocities  at 
infinite  distance  from  the  interface  and  with  the  y  axis 
directed  vertically  upward.  In  addition,  we  define  the 
capillary  number  k,  the  wave  speed  parameter  p  ,  and 
the  density  parameter  p  by  the  relations 

k  =  4tt27 /pi^X2  ,  (1) 

p  =  2irc*/g\  ,  (2) 

p=p2/p,  .  (3) 

Here,  p,  and  p2  are,  respectively,  the  densities  in  the 
lower  and  upper  fluids. 

On  the  interface,  the  Bernoulli  equation  and  the  pres¬ 
sure  jump  due  to  surface  tension  yield 

£<«i-«i)+(i-p*+£r  £-£u-p>.  (4) 

where  qx  and  q2  are  the  magnitudes  of  the  velocity  in  the  * 
lower  and  upper  fluids.  R  is  the  radius  of  curvature  of  .  f  A 
the  interface  counted  positive  when  the  center  of  curva¬ 
ture  lies  on  the  side  of  the  lower  fluid.  The  choice  of 
the  constant  in  Eq.  (4)  fixes  the  origin  of  y  as  the  undis¬ 
turbed  level  for  which  the  curvature  is  zero  and  qx  =q2 
=  1. 

We  shall  consider  z  =  x  +  iy  as  an  analytic  function  of 
the  complex  variables  fx  =  <t>\  +  Wi  and  /ts$s+j$a» 
respectively,  on  the  half  planes  $x  «  0  and  *  0.  Then, 
the  interface  can  be  represented  parametrically  in  two 
different  ways,  namely,  *($>,, 0.),  y(<t>\,0.)  and  x(4>2,0.), 
yU> which  we  shall  write  and  *,($>*), 

To  determine  the  shape  of  the  interface  we  note  that 
dx/d<px  +  idy/d<px  - 1  and  3x/3<>2  +t'd y/d<p2  - 1  vanish, 
respectively,  at  ilx  =  -  «  and  tf2  =  +«.  Therefore,  their 
real  parts  on  the  lines  ^  «0  and  </>2  =  0  are  the  Hilbert 
transforms  of  their  Imaginary  parts; 
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We  now  use  the  assumed  periodicity  and  symmetry  of 
the  surface  to  rewrite  (5)  and  (6)  in  the  forms 


(7) 

^-=1  +  j  [cotir(^j -4>2)“COtff(0j  +  <t>2j\d<Pi  • 

(8) 

Since  <px  *  <p 2  at  two  points  in  contact  on  the  interface, 
we  define  the  function 


,  (9) 

by  the  relations 

(<*>,)= *i*(<M  ,  do) 

ytM>t)=:vi*M>i)]  ,  (U) 

h(0)  =  0.  (12) 


Thus ,  we  have 

dx2  dxt  ( dh\~ 1 
d<P2  -d<px\d<fix)  * 

fdk  Y 

d<t> 2  d<t> j  \d<t>x  ) 


(13) 

(14) 


Substituting  Eqs.  (13)  and  (14)  into  Eq.  (8)  we  obtain 

C^r) "  - 1 + X’ ' ' 


-cot niM0;)  +  M0,)]}rf0(  •  (IS) 


Finally,  we  rewrite  condition  (4)  as 


JL 

4  n 


jc/dx^  .  rfjii  \[/£iY  +/2ztVT 
^  d*r /iWt)  w«/j 


a-i-i/2 


■fa'1-"1- 


u«i 


In  addition  to  the  parameters  k,  p,  andp,  a  given  wave 
is  characterized  by  a  third  parameter  which  is  a  mea¬ 
sure  of  the  wave  amplitude.  We  choose  this  parameter 
to  be  the  steepness  s  defined  by 

s=y(0)-y(ir).  (17) 


Dimensional  analysis  implies  that  a  functional  relation¬ 
ship  should  exist  among  the  four  numbers  *,p ,p,  and  s. 
We  fix  three  of  these  parameters  and  we  seek  three 
functions 

*,(♦,).  »nd*(<*>|).  el 0,1/21  , 

and  a  value  for  the  fourth  parameter  that  simultaneously 
satisfy  Eqs.  (7),  (15),  and  (16).  In  the  present  paper, 
we  fix  «,p,  and  s#  and  find  p  as  part  of  the  solution. 


ill.  NUMERICAL  ANALYSIS 


ence  method.  Introducing  a  uniform  mesh  we  have 

4>\={i -D/2N,  «  =  1 . N  + 1.  (18) 

Since  the  wave  is  symmetrical,  we  have 


Thus,  the  unknown  function  3y/9tf>  can  be  represented 
by  the  vector  of  dimension  N  - 1, 

Y'=(y,' , 

where 


We  now  define  the  midpoints  y4  =£(4>4  +  <£M),  t  =  l, 
N  and  we  represent  the  values  of 


dyt  dxx 


y(4> i)#  H<P{)t 


dh 

d<f>x  ' 


at  the  points  y,  by  the  vectors  yL  xm,  ym,  hmt  h*,  x" , 
and  yl. 


We  seek  to  satisfy  the  system  (7),  (15),  and  (16)  at 
the  points  y4.  The  integrals  in  Eqs.  (7)  and  (15)  are 
evaluated  at  the  points  y,  by  the  trapezoidal  rule  (which 
is  of  infinite  order  since  the  integrands  are  periodic). 
The  integration  is  over  the  points  4>{.  The  singularities 
of  the  Cauchy  principal  values  are  automatically  taken 
into  account  since  the  quadrature  is  symmetrical  with 
respect  to  the  singularities.  Thus,  we  obtain  from  Eqs. 
(7)  and  (15), 

x'..=  l+AY'  ,  (19) 

x'*=lu  +B(h)Y'  .  (20) 


Here,  A  is  a  known  matrix  and  B  a  matrix  whose  ele¬ 
ments  are  nonlinear  functions  of  the  vector 

*=[*(*?) . *(*?)]  • 

The  vector  y„  is  expressed  in  terms  of  Y'  by  a  sixth- 
order  quadrature  formula.  Thus, 

y*=Yo+CY\  (21) 

where  C  is  a  known  matrix.  The  elevation  Y0  of  the  in¬ 
terface  at  0,  =0  has  to  be  found  as  part  of  the  solution. 
Next,  we  express  xl,y«,h,  h'mt  and  y*  in  terms  of  x*, 
Y',  and  h„  by  sixth-order  interpolation  and  difference 
formulae. 


Substituting  Eqs.  (19)  and  (21)  into  Eq.  (15)  at  the 
points  yl#  i«l,  .  .  ,  ,  N  - 1  and  Eqs.  (16)  at  the  points 
y4 ,  i  *  1 ,  .  . .  ,N  we  obtain  a  system  of  2N  - 1  nonlinear 
algebraic  equations  for  the  2N  + 1  unknowns  Y' ,  hmtY ,, 
and  p .  Relations  (12)  and  (17)  provide  two  extra  equa¬ 
tions.  Thus,  tor  given  values  of  nc,p,  and  $  we  have 
a  system  of  2JV  +  1  equations  with  2^  +  1  unknowns.  This 
system  is  solved  by  Newton’s  Iterations. 


IV.  DISCUSSION  OF  RESULTS 
A.  Pure  gravity  waves 


We  seek  a  numerical  solution  of  the  integro-differen-  Before  proceeding  to  the  general  case  where  both 
tial  system  of  Eqs.  (7),  (15),  and  (16)  by  a  finite  differ-  gravity  and  surface  tension  are  taken  into  account,  we 
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shall  consider  the  limiting  case  o!  pure  gravity  waves. 

In  a  recent  paper,  Holyer1®  computed  accurate  solu¬ 
tions  for  pure  gravity  interfacial  waves  up  to  a  maxi¬ 
mum  value  of  the  steepness  at  which  the  profile  be¬ 
comes  vertical  at  some  point.  These  solutions  were 
obtained  by  calculating  the  coefficients  in  Stokes' 
expansion  on  a  computer  and  than  recasting  the  result¬ 
ing  high-order  polynomials  as  Padd  approximants. 

The  numerical  procedure  of  Sec,  in  was  used  to  com¬ 
pute  gravity  waves  with  p  =  0.1.  Preliminary  computa¬ 
tions  showed  that  an  accurate  solution  for  s  >  0.1  could 
not  be  computed  with  equal  increments  in  the  velocitv 
potential  4>t.  The  reason  is  that  the  spacing  of  the 
mesh  points  becomes  too  sparse  near  the  crest  where 
the  velocity  is  small.  A  similar  situation  was  pre¬ 
viously  encountered  by  Chen  and  Saffman, 8  Schwartz 
and  Vanden-Broeck,7  and  Vanden-Broeck  and  Schwartz8 
in  the  case  of  surface  waves.  These  authors  found  that 
this  difficulty  could  easily  be  overcome  by  concentrat¬ 
ing  the  mesh  points  near  the  crest  by  an  appropriate 
change  of  variable.  In  the  present  problem,  we  intro¬ 
duce  a  new  variable  0  by  the  relation 

J3=M4>i)  .  (22) 

Here,  the  functions  is  defined  by  Eqs.  (10)-(12).  The 
numerical  scheme  was  then  used  with  equal  increments 
in  the  new  variable  0.  In  Table  I,  we  present  values  of 
p  for  different  values  of  the  steepness  s  computed  with 
iV  =  15,  20,  and  25.  The  values  foriV  =  25  have  con¬ 
verged  to  five  decimal  places  for  s  <  0.13  and  to  three 
decimal  places  for  s  >  0.13.  These  values  agree  with 
those  presented  graphically  by  Holyer.10  Thus,  the 
validity  of  the  use  of  the  Pad£  approximant  method  as 
applied  to  gravity  interfacial  waves  is  confirmed.  A 
few  typical  profiles  are  shown  in  Fig.  1.  For  s~  0.222, 
the  profile  becomes  vertical  at  a  small  distance  from 
the  crest.  It  is  worthwhile  mentioning  that  our  num¬ 
erical  scheme  appears  to  be  more  efficient  than  the 
Padd  table  method  since  the  highest  wave  presented  by 
Holyer  corresponds  to  s  =  0.19. 

B.  Gravity-capillary  waves 

Some  insight  into  the  problem  can  be  gained  by  con¬ 
sidering  a  solution  in  the  form  of  a  Stokes'  expansion. 
Thus,  we  seek  a  solution  as  a  Fourier  expansion  in 
the  horizontal  coordinate  by  assuming  that  the  nth 
Fourier  coefficient  is  nth  order  in  the  amplitude.  In 
the  first  approximation,  the  waves  are  linear  sine  waves 


TABLE  I.  Values  ofp  for  0.05*  s*  0.22,  k  *0  and  p  =0.1. 


X 

N  =  l  S 

tf*20 

=  25 

0.05 

0.835097 

0.835108 

0.835107 

0.095493 

0.880605 

0.880620 

0.860624 

0.127324 

0.930  479 

0.930483 

0.930484 

0.159155 

0.996  561 

0.996  498 

0.996480 

0.190  986 

1.082  550 

1.082  075 

1.081  947 

0.206  90 

1.137 163 

1.135659 

1.135 199 

0.21 

1.149  422 

1.147491 

1.146  868 

0.22 

1.196  938 

1.192  921 

1.190  878 

o.i 

0. 


-0.1 

-0.2 

H - 1 - 1 - 1 - 1 - 1 - 

0  0.1  0.2  0.3  0.4  0.5 

FIG.  1,  Computed  wave  profiles  for  k  =  0,  p=0.1,  and  s  =  0.05p 
0.127,  0.2215. 


and  the  dispersion  relation  is  given  by 


M 


—  Lz£.  + _ L_ 

1+P  1-p* 


(23) 


It  can  easily  be  verified  by  substituting  Eqs.  (1)  and  (2) 
into  Eq.  (23)  that  linear  waves  of  wavelength  X  and  x/n 
travel  with  the  same  speed  c  if 


k  =  (1 -p)/n  , 


(24) 


where  n  is  an  integer  greater  than  1.  This  nonunique¬ 
ness  in  the  first  approximation  implies  that  some  of  the 
coefficients  in  the  Stokes'  expansion  will  become  infin¬ 
ite  when  k  assumes  one  of  the  values  (24).  In  the  par¬ 
ticular  case  p  =  0,  these  coefficients  have  been  com¬ 
puted  to  fifth  order  by  Wilton11  and  to  100th  order  by 
Hogan.5  Solutions  corresponding  to  the  critical  values 
(24)  can  be  found  by  revoking  Stokes'  hypothesis  and 
reordering  the  terms  of  the  expansion.  For  example, 
Wilton11  found  two  solutions  for  k  =  |  and  p  =  0,  i.e., 
for  the  critical  value  (24)  corresponding  to  «  =  2.  The 
numerical  work  of  Schwartz  and  Vanden-Broeck7  shows 
clearly  that  these  two  solutions  are,  in  fact,  members 


- 1 - 1 - 1 - 1 - 1 - * 

0.3  0.4  0.5  0.6  0,8 

FIG.  2.  Variation  of  speed  parameter  pwtth  capillary  number 
for  s  *  0.03  and  p  =  0.1. 
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-0.05 


-0.075' 


-0,1 25| 

- 1 - 1 - 1 - i - 1 - 1 - 

0  0.1  0.2  0.3  0.4  0.5 

FIG.  3.  Two  wave  profiles  for  *=0.45,  p=0.1,  and  s  =  0.03. 

of  two  different  families  labeled  1  and  2.  Families  1 
and  2  agree,  respectively,  with  the  Stokes*  expansion 
for  *  >  iand  i<  *<  5. 

In  the  present  section,  we  use  the  numerical  procedure 
of  Sec.  in,  to  compute  the  equivalent  of  these  two  fam¬ 
ilies  in  the  case  p  =  0.1.  The  critical  value  corres¬ 
ponding  to  w  =  2  in  Eq.  (24)  is  then  *  =  0.45.  We  started 
the  iterations  at  *  =  0.7  and  s  =0.03.  A  sine  wave  was 
used  as  an  initial  guess.  The  numerical  scheme  was 
found  to  converge  to  a  wave  of  family  1.  This  result 
could  be  expected  since  family  1  agrees  with  the  Stokes’ 
expansion  for  *  >  0.45.  This  solution  was  used  to  com¬ 
pute  the  solution  for  a  smaller  value  of  *  and  so  on. 
Family  2  was  computed  in  a  similar  way  by  starting 
the  iterations  at  *  =  0.4  with  a  sine  wave.  The  values 
of  the  speed  parameter  p  versus  *  for  s  =  0.03  are  pre¬ 
sented  in  Fig.  2.  Also  shown  is  the  infinitesimal  wave 
solution  (23).  The  two  families  cguld  be  extended  to 
smaller  and  larger  values  of  *.  This  was  not  done  since 
our  purpose  is  simply  to  describe  the  behavior  of  the 
solutions  near  *  =  0.45.  In  Fig.  3,  we  present  the  pro¬ 
files  of  the  waves  of  families  1  and  2  for  *  =  0.45  and 
s  =0.03.  They  are  the  direct  extensions  of  the  well- 
known  Wilton  ripples  for  interfacial  waves  with  p  =  0.1. 
In  Fig.  4,  we  present  a  steep  profile  of  a  wave  of  family 
I  for  *  =  0.45  and  s  =  0.15.  Although  indiscernible  in 
the  figure  the  wave  retains  a  slight  crest  dimple.  Near 
the  trough  the  slope  of  the  profile  varies  rapidly  and  is 
equal  to  zero  at  x  =  0,5.  We  were  unable  to  compute 
waves  of  higher  steepness  for  *  =  0.45.  The  numerical 


0  0.1  •  0.2  0.3  0.4  0.5 

—I - 1 - 1 - 1 - 1 - f- 

0  i 


-0.2 


FIG.  4.  Steep  profile  for  *  *  0.45. 

scheme  is  limited  by  the  small  values  of  the  velocity 
at  the  trough  in  the  upper  fluid.  We  expect  the  limiting 
profile  to  exhibit  a  small  trapped  bubble  at  the  trough 
with  a  stagnation  point  at  the  point  of  contact  in  the 
upper  fluid. 
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